#######################################################################
#Name of code file: mobilization_analysis.R

#Data In: st_panel.RDS

#Data Out: Figures 5, A10; Tables 2, A9, A10, A11, A12, A13, A14
######################################################################

# ----------------------------------------------------------------------
# load all relevant packages
# ----------------------------------------------------------------------

## If you don't have packages installed, uncomment and install
## install.packages("tidyverse")
## install.packages("ggplot2")
## install.packages("stargazer")
## install.packages("xtable")
## install.packages("texreg")
## install.packages("readstata13")
## install.packages("estimatr")
## install.packages("naniar")
## install.packages("lfe")
## install.packages("multiwayvcov")
## install.packages("lmtest")
## install.packages("xlsx")
## install.packages("readxl")
## install.packages("lubridate")
## install.packages("panelView")
## install.packages("zoo")
## Note that the Zelig package, which we use for making first differences plots,
## has been removed from CRAN - in order to install, you should first download
## the latest version from the Archive (5.1.7 as of this script) and manually
## install it using the following code - the first line should be the filepath
## to the .tar.gz file that you downloaded, the second line uses this filepath
## to install the package
## zeligpackage <- "~/Downloads/Zelig_5.1.7.tar.gz"
## install.packages(zeligpackage, repos = NULL, type = "source")
## install.packages("MASS")
## install.packages("sandwich")
## install.packages("lmtest")
## install.packages("bucky")
## install.packages("scales")

library("tidyverse")
library("ggplot2")
library("stargazer")
library("xtable")
library("texreg")
library("readstata13")
library("estimatr")
library("naniar")
library("lfe")
library("multiwayvcov")
library("lmtest")
# library("xlsx") Doesn't work
library("readxl")
library("lubridate")
library("panelView")
library("zoo")
library("Zelig")
library("MASS")
library("sandwich")
library("lmtest")
# library("bucky") 
library("scales")

## On the computer used to generate the results of this script, below is what is
## printed to the console when sessionInfo() is run
## sessionInfo()

## > sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.1

## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib

## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base

## other attached packages:
##  [1] scales_1.1.1       bucky_1.0.6        sandwich_3.0-1     MASS_7.3-54
##  [5] Zelig_5.1.7        survival_3.2-13    panelView_1.1.5    lubridate_1.8.0
##  [9] readxl_1.3.1       xlsx_0.6.5         lmtest_0.9-39      zoo_1.8-9
## [13] multiwayvcov_1.2.3 lfe_2.8-7.1        Matrix_1.3-4       naniar_0.6.1
## [17] estimatr_0.30.4    readstata13_0.10.0 texreg_1.37.5      xtable_1.8-4
## [21] stargazer_5.2.2    forcats_0.5.1      stringr_1.4.0      dplyr_1.0.7
## [25] purrr_0.3.4        readr_2.1.1        tidyr_1.1.4        tibble_3.1.6
## [29] ggplot2_3.3.5      tidyverse_1.3.1

## loaded via a namespace (and not attached):
##  [1] VGAM_1.1-5           colorspace_2.0-2     ellipsis_0.3.2
##  [4] class_7.3-19         visdat_0.5.3         fs_1.5.2
##  [7] rstudioapi_0.13      listenv_0.8.0        MatrixModels_0.5-0
## [10] prodlim_2019.11.13   fansi_0.5.0          xml2_1.3.3
## [13] codetools_0.2-18     splines_4.1.2        Formula_1.2-4
## [16] jsonlite_1.7.2       pROC_1.18.0          mcmc_0.9-7
## [19] caret_6.0-90         rJava_1.0-6          broom_0.7.10
## [22] dbplyr_2.1.1         geepack_1.3-2        compiler_4.1.2
## [25] httr_1.4.2           backports_1.4.0      assertthat_0.2.1
## [28] survey_4.1-1         cli_3.1.0            quantreg_5.86
## [31] tools_4.1.2          coda_0.19-4          gtable_0.3.0
## [34] glue_1.5.1           reshape2_1.4.4       Rcpp_1.0.7
## [37] carData_3.0-4        cellranger_1.1.0     vctrs_0.3.8
## [40] Amelia_1.8.0         nlme_3.1-153         conquer_1.2.1
## [43] iterators_1.0.13     timeDate_3043.102    gower_0.2.2
## [46] globals_0.14.0       xlsxjars_0.6.1       rvest_1.0.2
## [49] lifecycle_1.0.1      future_1.23.0        ipred_0.9-12
## [52] miscTools_0.6-26     hms_1.1.1            parallel_4.1.2
## [55] SparseM_1.81         gridExtra_2.3        rpart_4.1-15
## [58] stringi_1.7.6        foreach_1.5.1        AER_1.2-9
## [61] boot_1.3-28          lava_1.6.10          matrixStats_0.61.0
## [64] rlang_0.4.12         pkgconfig_2.0.3      lattice_0.20-45
## [67] recipes_0.1.17       tidyselect_1.1.1     parallelly_1.29.0
## [70] plyr_1.8.6           magrittr_2.0.1       R6_2.5.1
## [73] generics_0.1.1       DBI_1.1.1            pillar_1.6.4
## [76] haven_2.4.3          foreign_0.8-81       withr_2.4.3
## [79] nnet_7.3-16          abind_1.4-5          future.apply_1.8.1
## [82] modelr_0.1.8         crayon_1.4.2         car_3.0-12
## [85] utf8_1.2.2           tzdb_0.2.0           maxLik_1.5-2
## [88] grid_4.1.2           data.table_1.14.2    ModelMetrics_1.2.2.2
## [91] reprex_2.0.1         digest_0.6.29        MCMCpack_1.6-0
## [94] stats4_4.1.2         MatchIt_4.3.2        munsell_0.5.0
## [97] mitools_2.4

# ----------------------------------------------------------------------
# load ST clean data
# ----------------------------------------------------------------------

## Uncomment and set working directory to replication folder
## setwd("folder")
rm(list = ls())
setwd("C:/Users/User/Desktop/EUI/2nd term/Replicating Research/10/BSS replication files/Weiss et al original files/doc_replication")


st_final <- readRDS("data/st_final.RDS")

# ----------------------------------------------------------------------
# Merge ST data with Locality panel
# ----------------------------------------------------------------------

# replace NA with 0
final_data <- st_final %>%
    mutate(.,
           daily_join = ifelse(is.na(daily_join), 0, daily_join),
           join_binary = ifelse(daily_join >0,1,0),
           week = lubridate::week(ymd(st_final$date)),
           year = lubridate::year(ymd(st_final$date)),
           month = lubridate::month(ymd(st_final$date))) %>%
    mutate(year_fac = factor(year),
           month_fac = factor(month),
           week_fac = factor(week),
           lcode_fac = factor(lcode))

# ----------------------------------------------------------------------
# Assign post deal date
# ----------------------------------------------------------------------

final_data <- final_data %>%
 mutate(.,
        post = ifelse(date > "2020-01-27",1,0))

# ----------------------------------------------------------------------
# Mixed and Arab only data
# ----------------------------------------------------------------------

# Create subset of non-jewish localities

final_data_nj <- final_data %>%
 filter(.,
        non_jewish == 1)

sum(final_data_nj$daily_join, na.rm=T)
nrow(final_data_nj)

# ----------------------------------------------------------------------
# Create descriptive stats table
# ----------------------------------------------------------------------

disc_table_nj <- dplyr::select(final_data_nj,
                        c(daily_join, join_binary,
                          triangle, ext_tri, pop_2018, age0_19sef_pcnt_08,
                          age65sef_pcnt_08, age85sef_pcnt_08,
                          AcadmCert_pcnt_08, Worked_08, HousingDens_avg_08,
                          Vehicle1up_pcnt_08, ChldBorn_avg_08))
stargazer(as.data.frame(disc_table_nj),
          covariate.labels = c("Daily Join (Count)", "Daily Join (Binary)",
                               "Triangle", "Extended Triangle",  "Population 2018",
                               "Perc. Age 0-19", "Perc. Age 65+", "Perc. Age 85+",
                               "Perc. Academic", "Perc. Employed", "Housing Density",
                               "HH with Vehicle", "Average Children per Woman"),
          title = "Descriptive Statistics - Non Jewish Localities",
          label = "tab:dsc_st",
          style = "qje",
          omit.summary.stat = c("p25", "p75"),
          notes = c("All variables following 'Perc. Age 0-19' are from the 2008 census."),
          notes.append = T,
          notes.align = "l",
          out="tables/table_A9.tex")

# ----------------------------------------------------------------------
# Basic DiD
# ----------------------------------------------------------------------

## Table 2 model 1: main specification
m1 <- felm(join_binary ~ triangle * post + pop_2018 | 0 | 0 | lcode,
           data = final_data_nj,
           exactDOF = TRUE,
           keepCX = TRUE)
summary(m1)

## Table 2 model 2: main specification with time and locality fixed effects
m2 <- felm(join_binary ~ triangle * post | lcode + year + month + week | 0 | lcode,
           data = final_data_nj,
           exactDOF = TRUE,
           keepCX = TRUE)
summary(m2)

## Table 2 model 3: alternative treatment with time and locality fixed effects
m3 <- final_data_nj %>%
    mutate(.,
           triangle = ext_tri) %>%
    felm(join_binary ~ triangle * post | lcode + year + month + week | 0 | lcode,
         data = .,
         exactDOF = TRUE,
         keepCX = TRUE)
summary(m3)

## Table 2 model 4: full sample with time and locality fixed effects
m4 <- felm(join_binary ~ triangle * post | lcode + year + month + week | 0 | lcode,
           data = final_data,
           exactDOF = TRUE,
           keepCX = TRUE)
summary(m4)

## Remove SEs that are zero when coefficient is NA
se_list <- lapply(list(m1, m2, m3, m4), function(x) coefficients(summary(x))[, "Cluster s.e."])
se_list[[2]][1] <- NA
se_list[[3]][1] <- NA
se_list[[4]][1] <- NA

## Table for paper
stargazer(m1, m2, m3, m4,
          out= "tables/table2.tex",  style="qje",
          title="Deal of the Century Effect on Mobilization",
          keep = c("triangle", "post", "triangle:post"),
          dep.var.labels = c("Mobilization"),
          covariate.labels = c("Triangle", "Post", "Triangle * Post"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          se = se_list,
          omit.table.layout = "n",
          column.sep.width = "-6pt",
          star.cutoffs = NULL,
          label = "tab:main_did_st",
          star.char = c("", "", ""),
          add.lines = list(
              c("Population Control", "Yes", "No", "No", "No"),
              c("Week FE", "No", "Yes", "Yes", "Yes"),
              c("Month FE", "No", "Yes", "Yes", "Yes"),
              c("Year FE", "No", "Yes", "Yes", "Yes"),
              c("Locality FE", "No", "Yes", "Yes", "Yes"),
              c("Cluster",  "Locality", "Locality", "Locality", "Locality"),
              c("Sample",  "Non-Jewish", "Non-Jewish", "Non-Jewish", "Full"),
              c("Treatment",  "10 Localities", "10 Localities", "16 Localities",
                "10 Localities"),
              c("Pre-Register", "No", "No", "No", "No")))

## Analysis for appendix with combined locality-year-month-week effects
modnj_lfe_x_time <- felm(join_binary ~ triangle * post | paste(lcode, year, month, week) | 0 | lcode,
                         data = final_data_nj,
                         exactDOF = TRUE,
                         keepCX = TRUE)
summary(modnj_lfe_x_time)
modfull_lfe_x_time <- felm(join_binary ~ triangle * post | paste(lcode, year, month, week) | 0 | lcode,
                           data = final_data,
                           exactDOF = TRUE,
                           keepCX = TRUE)
summary(modfull_lfe_x_time)

## Table for paper, edit manually to remove extraneous 0s for SEs
stargazer(modnj_lfe_x_time, modfull_lfe_x_time,
          out= "tables/table_A10.tex",  style="qje",
          title="Deal of the Century Effect on Mobilization",
          keep = c("post", "triangle:post"),
          dep.var.labels = c("Mobilization"),
          covariate.labels = c("Post", "Triangle * Post"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          omit.table.layout = "n",
          column.sep.width = "-6pt",
          star.cutoffs = NULL,
          label = "tab:st_lfextime",
          star.char = c("", "", ""),
          add.lines = list(c("Locality-Week Fixed Effects", "Yes", "Yes"),
                           c("Sample", "Non-Jewish", "Full")))

# ----------------------------------------------------------------------
# Negative binomial models
# ----------------------------------------------------------------------

## Negative binomial models
nb1 <- glm.nb(daily_join ~ triangle * post + pop_2018, data = final_data_nj)
summary(nb1)
nb2 <- glm.nb(daily_join ~ triangle * post + pop_2018, data = final_data)
summary(nb2)

## Negative binomial model but in Zelig so that I can easily plot first
## differences
set.seed(84058)
nb1_alt <- zelig(daily_join ~ triangle * post + pop_2018, model = "negbin",
                 data = final_data_nj)
x0 <- setx(nb1_alt, triangle = 1, post = 0)
x1 <- setx(nb1_alt, triangle = 1, post = 1)
sim_out <- sim(nb1_alt, x = x0, x1 = x1)
summary(sim_out)
pdf("plots/figure_A10.pdf", width = 5.75, height = 7)
plot(sim_out)
dev.off()

## Output in table
stargazer(nb1, nb2,
          out= "tables/table_A14.tex",  style="qje",
          title="Deal of the Century Effect on Mobilization, Negative Binomial Models",
          keep = c("triangle", "post", "triangle:post"),
          dep.var.labels = c("Mobilization"),
          covariate.labels = c("Triangle", "Post", "Triangle * Post"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          omit.table.layout = "n",
          column.sep.width = "-6pt",
          star.cutoffs = NULL,
          label = "tab:negbin",
          star.char = c("", "", ""),
          add.lines = list(c("Population Controls", "Yes", "Yes"),
                           c("Sample", "Non-Jewish", "Full")))

# ----------------------------------------------------------------------
# Analysis with restricted date range
# ----------------------------------------------------------------------

m1altd <- felm(join_binary ~ triangle*post |0| 0 |lcode,
               data = final_data_nj,
               subset = date > "2019-01-01",
               exactDOF=TRUE,
               keepCX = TRUE)
m2altd <- felm(join_binary ~ triangle*post |week + month + year| 0 |lcode,
               data = final_data_nj,
               subset = date > "2019-01-01",
               exactDOF=TRUE,
               keepCX = TRUE)
m3altd <- felm(join_binary ~ triangle*post + pop_2018|week + month + year| 0 |lcode,
               data = final_data_nj,
               subset = date > "2019-01-01",
               exactDOF=TRUE,
               keepCX = TRUE)
m4altd <- final_data_nj %>%
    mutate(.,
           triangle = ext_tri) %>%
    felm(join_binary ~ triangle*post + pop_2018
         |week + month + year| 0 |lcode,
         data = .,
         subset = date > "2019-01-01",
         exactDOF=TRUE,
         keepCX = TRUE)
m5altd <- felm(join_binary ~ triangle*post |week + month + year| 0 |lcode,
               data = final_data,
               subset = date > "2019-01-01",
               exactDOF=TRUE,
               keepCX = TRUE)
stargazer(m1altd, m2altd, m3altd, m4altd, m5altd,
          out= "tables/table_A13.tex",  style="qje",
          title="Deal of the Century Effect on Mobilization, 2019--2020",
          keep = c("triangle", "post", "triangle:post"),
                                        #omit = c("Constant", "cycleThird", "cycleSecond", "pop_2018", "AcadmCert_pcnt_08"),
          dep.var.labels = c("Mobilization"),
          covariate.labels = c("Triangle", "Post", "Triangle*Post"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          omit.table.layout = "n",
          column.sep.width = "-6pt",
          star.cutoffs = NULL,
          label = "tab:main_did_st_altd",
          star.char = c("", "", ""),
          add.lines = list(
              c("Week FE",  "No", "Yes", "Yes", "Yes", "Yes"),
              c("Month FE",  "No", "Yes", "Yes", "Yes", "Yes"),
              c("Year FE",  "No", "Yes", "Yes", "Yes", "Yes"),
              c("Pop Control",  "No", "No", "Yes","Yes", "No"),
              c("Cluster",  "Locality", "Locality", "Locality", "Locality","Locality"),
              c("Sample",  "Non-Jewish", "Non-Jewish", "Non-Jewish", "Non-Jewish",
                "Full"),
              c("Treatment",  "10 Localities", "10 Localities", "10 Localities",
                "16 Localities", "10 Localities"),
              c("Pre-Register",  "No",  "No","No", "No", "No")))

# ----------------------------------------------------------------------
# With census covariates
# ----------------------------------------------------------------------

rb1 <- felm(join_binary ~ triangle*post +
                      AcadmCert_pcnt_08 + Worked_08
                     | week + month +year | 0 |lcode,
                     data = final_data_nj,
                     exactDOF=TRUE,
                     keepCX = TRUE)
summary(rb1)

rb2 <- felm(join_binary ~ triangle*post  +
                      AcadmCert_pcnt_08 + Worked_08 +
                      HousingDens_avg_08 + Vehicle1up_pcnt_08 + ChldBorn_avg_08
                     | week + month +year | 0 |lcode,
                     data = final_data_nj,
                     exactDOF=TRUE,
                     keepCX = TRUE)
summary(rb2)

rb3 <- felm(join_binary ~ triangle*post  +
                      AcadmCert_pcnt_08 + Worked_08 +
                      HousingDens_avg_08 + Vehicle1up_pcnt_08 + ChldBorn_avg_08+
                      age0_19sef_pcnt_08 + age65sef_pcnt_08 + age85sef_pcnt_08
                     | week + month +year | 0 |lcode,
                     data = final_data_nj,
                     exactDOF=TRUE,
                     keepCX = TRUE)
summary(rb3)

rb4 <- felm(join_binary ~ triangle*post +
             AcadmCert_pcnt_08 + Worked_08
            | week + month +year | 0 |lcode,
            data = final_data,
            exactDOF=TRUE,
            keepCX = TRUE)
summary(rb4)

rb5 <- felm(join_binary ~ triangle*post  +
             AcadmCert_pcnt_08 + Worked_08 +
             HousingDens_avg_08 + Vehicle1up_pcnt_08 + ChldBorn_avg_08
            | week + month +year | 0 |lcode,
            data = final_data,
            exactDOF=TRUE,
            keepCX = TRUE)
summary(rb5)

rb6 <- felm(join_binary ~ triangle*post  +
             AcadmCert_pcnt_08 + Worked_08 +
             HousingDens_avg_08 + Vehicle1up_pcnt_08 + ChldBorn_avg_08+
             age0_19sef_pcnt_08 + age65sef_pcnt_08 + age85sef_pcnt_08
            | week + month +year | 0 |lcode,
            data = final_data,
            exactDOF=TRUE,
            keepCX = TRUE)
summary(rb6)

rb3add <- felm(join_binary ~ triangle * post
            | week + month + year_fac + AcadmCert_pcnt_08:year_fac +
              Worked_08:year_fac + HousingDens_avg_08:year_fac +
              Vehicle1up_pcnt_08:year_fac + ChldBorn_avg_08:year_fac +
              age0_19sef_pcnt_08:year_fac + age65sef_pcnt_08:year_fac +
              age85sef_pcnt_08:year_fac
            | 0 | lcode,
            data = final_data_nj,
            exactDOF = TRUE,
            keepCX = TRUE)
summary(rb3add)

rb6add <- felm(join_binary ~ triangle * post
            | week + month + year_fac + AcadmCert_pcnt_08:year_fac +
              Worked_08:year_fac + HousingDens_avg_08:year_fac +
              Vehicle1up_pcnt_08:year_fac + ChldBorn_avg_08:year_fac +
              age0_19sef_pcnt_08:year_fac + age65sef_pcnt_08:year_fac +
              age85sef_pcnt_08:year_fac
            | 0 | lcode,
            data = final_data,
            exactDOF = TRUE,
            keepCX = TRUE)
summary(rb6add)

## Output in table
stargazer(rb1, rb2, rb3, rb3add, rb4, rb5, rb6, rb6add,
          out= "tables/table_A11.tex",  style="qje",
          title="Deal of the Century Effect on Mobilization (2008 Census Covariates)",
          order = c(1,2,11,3:20),
          omit = c("Constant"),
          font.size = "scriptsize",
          dep.var.labels = c("Mobilization"),
          covariate.labels = c("Triangle", "Post", "Triangle*Post",
                               "Perc. Academic", "Perc. Employed", "Housing Density",
                               "HH with Vehicle", "Average Children per Women",
                               "Perc. Age 0-19", "Perc. Age 65+",
                               "Perc. Age 85+"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
        #  font.size = "tiny",
          omit.table.layout = "n",
          column.sep.width = "-6pt",
          star.cutoffs = NULL,
          label = "tab:rob_cov_st",
          star.char = c("", "", ""),
          add.lines = list(
              c("Week FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
              c("Month FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
              c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
              c("Year * Census FE", "No", "No", "No", "Yes", "No", "No", "No", "Yes"),
              c("Cluster", "Locality", "Locality", "Locality", "Locality",
                "Locality", "Locality", "Locality", "Locality"),
              c("Sample", "Non-Jewish", "Non-Jewish", "Non-Jewish", "Non-Jewish",
                "All", "All", "All", "All"),
              c("Pre-Register", "No", "No", "No", "No", "No", "No", "No", "No")))

# ----------------------------------------------------------------------
# Logit
# ----------------------------------------------------------------------

logit_1 <- glm(join_binary ~ triangle*post,
              family = "binomial",
              data = final_data_nj)
logit_1_cluster <-robustify(logit_1, cluster= lcode)

logit_2 <- glm(join_binary ~ triangle*post +
                as.factor(week) + as.factor(month) + as.factor(year),
               family = "binomial",
               data = final_data_nj)
logit_2_cluster <-robustify(logit_2, cluster= lcode)


logit_3 <- glm(join_binary ~ triangle*post + pop_2018 +
                as.factor(week) + as.factor(month) + as.factor(year),
               family = "binomial",
               data = final_data_nj)
logit_3_cluster <-robustify(logit_3, cluster= lcode)

## Output in table
stargazer(logit_1_cluster, logit_2_cluster, logit_3_cluster,
          out= "tables/table_A12.tex",  style="qje",
          title="Deal of the Century Effect on Mobilization (Logit)",
          keep = c("triangle", "post", "triangle:post"),
          dep.var.labels = c("Mobilization"),
          covariate.labels = c("Triangle", "Post", "Triangle*Post"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq", "aic", "ll"),
          ci=FALSE,
          omit.table.layout = "n",
          column.sep.width = "-6pt",
          star.cutoffs = NULL,
          label = "tab:rob_st_logit",
          star.char = c("", "", ""),
          add.lines = list(
              c("Week FE",  "No", "Yes", "Yes"),
              c("Month FE",  "No", "Yes", "Yes"),
              c("Year FE",  "No", "Yes", "Yes"),
              c("Pop Control",  "No", "No", "Yes"),
              c("Cluster",  "Locality", "Locality", "Locality"),
              c("Sample",  "Non-Jewish", "Non-Jewish", "Non-Jewish"),
              c("Pre-Register",  "No", "No", "No")))

# ----------------------------------------------------------------------
# Parallel trends
# ----------------------------------------------------------------------

# Group ST data by locality and weekid
pt <- group_by(final_data_nj, triangle, year, month)

# Collapse ST data by locality and date
sign_up <- summarize(pt,
                     daily_join = sum(daily_join, na.rm = T))
sign_up <- as.data.frame(sign_up)
sign_up$date <- as.yearmon(paste(sign_up$year,sign_up$month), "%Y %m")
sign_up <- mutate(sign_up,
                  triangle = ifelse(triangle==1, "Triangle", "Non-Triangle"))
sign_up$triangle <- factor(sign_up$triangle, c("Triangle", "Non-Triangle"))

## Plot
ggplot(sign_up, aes(x = date, y = daily_join, color = as.factor(triangle),
                    shape = as.factor(triangle))) +
    geom_line(aes(color=as.factor(triangle)))+
    geom_point(aes(color=as.factor(triangle))) +
    theme(axis.text.x=element_text(size=10))+
    geom_vline(xintercept=as.yearmon("2020-01-28"), linetype="dotdash",
               color = "grey", size=0.6) +
    annotate("text", x = as.yearmon("2019-01-28"), y = 1400, label = "PRE",
             family = "Times", size = 3) +
    annotate("text", x = as.yearmon("2020-04-01"), y = 1400, label = "POST",
             family = "Times", size = 3) +
    scale_color_manual(values = c("firebrick3", "dodgerblue3")) +
    scale_shape_manual(values = c("triangle", "circle")) +
    scale_y_continuous(labels = label_comma()) +
    labs(color = "", shape = "",
         x = "",
         y = "Monthly Mobilization")+
    theme(text = element_text(size = 12, family = "Times"),
          legend.key=element_blank(),
          panel.grid.major = element_blank(),
          axis.text.x = element_text(size = 12),
          plot.caption = element_text(size = 12, family = "Times",hjust = -.02),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(colour = "black"))
ggsave("plots/figure5.pdf", width = 6.5, height = 4)





#################### ---------- BSS -----------------------###############################################################################
###
###
###
### ---- Mobilization

library(fixest)

## Event study
# days
df_es <- final_data %>%
  mutate(
    event_time_day = as.numeric(date - as.Date("2020-01-28"))
  ) %>%
  filter(event_time_day >= -60 & event_time_day <= 60)  # Keep only ±60 days
event_model_day <- feols(
  join_binary ~ i(event_time_day, triangle, ref = -1) | lcode + event_time_day,
  cluster = "lcode",
  data = df_es
)

png("plots/BSS_figure_3.png", width = 6.5, height = 4, units = "in", res = 300)
iplot(event_model_day,
      xlab = "Days since treatment",
      ylab = "Effect on Joining Movement",
      main = "Event Study: ±60 Days Around Treatment",
      col = "firebrick3")
dev.off()

# months
df <- final_data_nj %>%
  arrange(year, month) %>%
  mutate(
    month_seq = (year - 2020) * 12 + (month - 1)
  )

did_3 <- feols(
  join_binary ~ i(month_seq, triangle, ref = -1) | lcode + month_seq, 
  cluster = "lcode", 
  data = df
)

png("plots/BSS_figure_4.png", width = 6.5, height = 4, units = "in", res = 300)
iplot(did_3,
      xlab = "Months since treatment (Jan 2020)",
      ylab = "Effect on Joining Movement",
      main = "Event Study: Monthly Time Scale Centered on Treatment",
      col = "firebrick3")
dev.off()



## TABLE 6 ---------------------
# Model 1
# Their model m2: 10 treated localities with time and locality fixed effects

# Model 2
# Their model m3: 16 treated localities with time and locality fixed effects

# Model 3
# Our model m5: 6 (not-mentioned) treated localities, 10 mentioned ones included as controls
m5 <- final_data_nj %>%
  mutate(triangle = ifelse(ext_tri == 1 & triangle == 0, 1, 0)) %>%
  felm(join_binary ~ triangle * post | lcode + year + month + week | 0 | lcode,
       data = ., exactDOF = TRUE, keepCX = TRUE)
summary(m5)

# Model 3
# Our model m6: 6 (not-mentioned) treated localities, 10 mentioned ones excluded from controls
m6 <- final_data_nj %>%
  mutate(triangle = ifelse(ext_tri == 1 & triangle == 0, 1, 0)) %>%
  filter(triangle == 1 | ext_tri == 0) %>%
  felm(join_binary ~ triangle * post | lcode + year + month + week | 0 | lcode,
       data = ., exactDOF = TRUE, keepCX = TRUE)
summary(m6)

## TABLE OUTPUT ##
models_mobilisation <- list(m2, m3, m5, m6)

stargazer(models_mobilisation,
          out= "tables/table6_BSS.tex",  style="qje",
          title="Deal of the Century Effect on Mobilization",
          keep = c("triangle", "post", "triangle:post"),
          dep.var.labels = c("Mobilization"),
          covariate.labels = c("Triangle", "Post", "Triangle * Post"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci=FALSE,
          se = se_list,
          omit.table.layout = "n",
          column.sep.width = "12pt",
          label = "tab:main_did_mobilisation",
          add.lines = list(
            c("Population Control", "No", "No", "No", "No"),
            c("Week FE", "Yes", "Yes", "Yes", "Yes"),
            c("Month FE", "Yes", "Yes", "Yes", "Yes"),
            c("Year FE", "Yes", "Yes", "Yes", "Yes"),
            c("Locality FE", "Yes", "Yes", "Yes", "Yes"),
            c("Cluster", "Locality", "Locality", "Locality", "Locality"),
            c("Sample",  "Non-Jewish", "Non-Jewish", "Non-Jewish", "No 10"),
            c("Treatment", "10 Loc.", "16 Loc.", "6 Loc.", "6 Loc."),
            c("Pre-Register", "No", "No", "No", "No")))






## Figure 5 
#visualising parallel trends, recoding 

library(lubridate)
library(scales)
library(dplyr)


final_data_nj$date <- as.Date(final_data_nj$date)

# after mid-2019, by week
pt <- final_data_nj %>%
  filter(date > as.Date("2019-07-01")) %>%  
  mutate(week = floor_date(date, "week")) %>%  
  group_by(triangle, week) %>%
  dplyr::summarize(daily_join = sum(daily_join, na.rm = TRUE))


pt <- mutate(pt,
             triangle = ifelse(triangle == 1, "Triangle", "Non-Triangle"))
pt$triangle <- factor(pt$triangle, levels = c("Triangle", "Non-Triangle"))

ggplot(pt, aes(x = week, y = daily_join, color = as.factor(triangle),
               shape = as.factor(triangle))) +
  geom_line(aes(color = as.factor(triangle))) +
  geom_point(aes(color = as.factor(triangle))) +
  theme(axis.text.x = element_text(size = 10)) +
  geom_vline(xintercept = as.Date("2020-01-28"), linetype = "dotdash",
             color = "grey", size = 0.6) +
  annotate("text", x = as.Date("2019-09-01"), y = 1400, label = "PRE",
           family = "Times", size = 3) +
  annotate("text", x = as.Date("2020-04-01"), y = 1400, label = "POST",
           family = "Times", size = 3) +
  scale_color_manual(values = c("firebrick3", "dodgerblue3")) +
  scale_shape_manual(values = c("triangle", "circle")) +
  scale_y_continuous(labels = label_comma()) +
  labs(color = "", shape = "",
       x = "",
       y = "Weekly Mobilization") +
  theme(text = element_text(size = 12, family = "Times"),
        legend.key = element_blank(),
        panel.grid.major = element_blank(),
        axis.text.x = element_text(size = 12),
        plot.caption = element_text(size = 12, family = "Times", hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))

ggsave("plots/BSS_figure_5.png", width = 6.5, height = 4)






# rerunning models with different operationalisation of the outcome variable

final_data_nj$join_binary_3plus <- ifelse(final_data_nj$daily_join >= 3, 1, 0)
final_data$join_binary_3plus <- ifelse(final_data$daily_join >= 3, 1, 0)

m1_thresh3 <- felm(join_binary_3plus ~ triangle * post + pop_2018 | 0 | 0 | lcode,
                   data = final_data_nj,
                   exactDOF = TRUE,
                   keepCX = TRUE)
summary(m1_thresh3)

## Table 2 model 2: main specification with time and locality fixed effects
m2 <- felm(join_binary_3plus ~ triangle * post | lcode + year + month + week | 0 | lcode,
           data = final_data_nj,
           exactDOF = TRUE,
           keepCX = TRUE)
summary(m2)

## Table 2 model 3: alternative treatment with time and locality fixed effects
m3 <- final_data_nj %>%
  mutate(.,
         triangle = ext_tri) %>%
  felm(join_binary_3plus ~ triangle * post | lcode + year + month + week | 0 | lcode,
       data = .,
       exactDOF = TRUE,
       keepCX = TRUE)
summary(m3)

## Table 2 model 4: full sample with time and locality fixed effects
m4 <- felm(join_binary_3plus ~ triangle * post | lcode + year + month + week | 0 | lcode,
           data = final_data,
           exactDOF = TRUE,
           keepCX = TRUE)
summary(m4)


# OUTPUT 
se_list <- list(
  summary(m1_thresh3)$coef[, "Cluster s.e."],
  summary(m2)$coef[, "Cluster s.e."],
  summary(m3)$coef[, "Cluster s.e."],
  summary(m4)$coef[, "Cluster s.e."]
)

se_list <- lapply(se_list, function(x) { x[x < 1e-6] <- NA; x })


stargazer(m1_thresh3, m2, m3, m4,
          out = "tables/table7_BSS.tex",
          style = "qje",
          title = "Deal of the Century Effect on Mobilization (3+ Subscriptions Threshold)",
          keep = c("triangle", "post", "triangle:post"),
          dep.var.labels = c("Mobilization (3+ joins)"),
          covariate.labels = c("Triangle", "Post", "Triangle * Post"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"),
          ci = FALSE,
          se = se_list,
          omit.table.layout = "n",
          column.sep.width = "-6pt",
          star.cutoffs = NULL,
          star.char = c("", "", ""),
          label = "tab:main_did_st_thresh3",
          add.lines = list(
            c("Population Control", "Yes", "No", "No", "No"),
            c("Week FE", "No", "Yes", "Yes", "Yes"),
            c("Month FE", "No", "Yes", "Yes", "Yes"),
            c("Year FE", "No", "Yes", "Yes", "Yes"),
            c("Locality FE", "No", "Yes", "Yes", "Yes"),
            c("Cluster", "Locality", "Locality", "Locality", "Locality"),
            c("Sample", "Non-Jewish", "Non-Jewish", "Non-Jewish", "Full"),
            c("Treatment", "10 Localities", "10 Localities", "16 Localities", "10 Localities"),
            c("Pre-Register", "No", "No", "No", "No")
          )
)







